## onefile=FALSE
## pdf(file = ifelse(onefile, "Rplots.pdf", "Rplot%03d.pdf"),
##         encoding="ISOLatin1",
##         width = 10, height = 8, onefile = TRUE, family = "Helvetica",
##         title = "R Graphics Output", fonts = NULL, version = "1.1")

##
## setting up plotting space
##
par(mfrow = c(1, 1))

##
## reading data and setting parameters for the nocascade (left) panel
##
data_1000_cascadeD0001 <- read.csv(file="../figure7_1000_cascadeD0001.csv",header=FALSE)
data_1000_cascadeD001 <- read.csv(file="../figure7_1000_cascadeD001.csv",header=FALSE)
data_1000_cascadeD01 <- read.csv(file="../figure7_1000_cascadeD01.csv",header=FALSE)
data_1000_cascadeD1 <- read.csv(file="../figure7_1000_cascadeD1.csv",header=FALSE)

data_1000_cascadeD0001 <- read.csv(file="../figure7_1000_cascadeD0001.csv",header=FALSE)
data_1000_cascadeD001 <- read.csv(file="../figure7_1000_cascadeD001.csv",header=FALSE)
data_1000_cascadeD01 <- read.csv(file="../figure7_1000_cascadeD01.csv",header=FALSE)
data_1000_cascadeD1 <- read.csv(file="../figure7_1000_cascadeD1.csv",header=FALSE)

data_2000_cascadeD0001 <- read.csv(file="../figure7_2000_cascadeD0001.csv",header=FALSE)
data_2000_cascadeD001 <- read.csv(file="../figure7_2000_cascadeD001.csv",header=FALSE)
data_2000_cascadeD01 <- read.csv(file="../figure7_2000_cascadeD01.csv",header=FALSE)
data_2000_cascadeD1 <- read.csv(file="../figure7_2000_cascadeD1.csv",header=FALSE)

data_2000_cascadeD0001 <- read.csv(file="../figure7_2000_cascadeD0001.csv",header=FALSE)
data_2000_cascadeD001 <- read.csv(file="../figure7_2000_cascadeD001.csv",header=FALSE)
data_2000_cascadeD01 <- read.csv(file="../figure7_2000_cascadeD01.csv",header=FALSE)
data_2000_cascadeD1 <- read.csv(file="../figure7_2000_cascadeD1.csv",header=FALSE)

minX.pow <- 3
maxX.pow <- 6
minY.pow <- -4
maxY.pow <- 0

xlab=expression(paste(n[0](cm)))

par(mar = c(6, 6, 3, 3))
plot(NULL, xlim=c(10^minX.pow, 10^maxX.pow), ylim=c(10^minY.pow, 10^maxY.pow),
     log = "xy", xaxt = "n",  yaxt = "n", type = "l", xlab=expression(paste(n[0](cm^{-3}))), 
     ylab=expression(paste(HDO/H[2], O)), col=1)

data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_1000_cascadeD0001[1,2:5]/data_1000_cascadeD0001[2,2:5]) )
lines(data1, col=1, lty=1)
points(data1, pch=8, cex=2, col=1, lty=1)
data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_1000_cascadeD001[1,2:5]/data_1000_cascadeD001[2,2:5]) )
lines(data1, col=1, lty=1)
points(data1, pch=5, cex=2, col=1, lty=1)
data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_1000_cascadeD01[1,2:5]/data_1000_cascadeD01[2,2:5]) )
lines(data1, col=1, lty=1)
points(data1, pch=2, cex=2, col=1, lty=1)
data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_1000_cascadeD1[1,2:5]/data_1000_cascadeD1[2,2:5]) )
lines(data1, col=1, lty=1)
points(data1, pch=22, cex=2, col=1, lty=1)

data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_2000_cascadeD0001[1,2:5]/data_2000_cascadeD0001[2,2:5]) )
lines(data1, col=2, lty=2)
points(data1, pch=8, cex=2, col=2, lty=1)
data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_2000_cascadeD001[1,2:5]/data_2000_cascadeD001[2,2:5]) )
lines(data1, col=2, lty=2)
points(data1, pch=5, cex=2, col=2, lty=1)
data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_2000_cascadeD01[1,2:5]/data_2000_cascadeD01[2,2:5]) )
lines(data1, col=2, lty=2)
points(data1, pch=2, cex=2, col=2, lty=1)
data1 <- t( rbind(c(1E3,1E4,1E5,1E6), data_2000_cascadeD1[1,2:5]/data_2000_cascadeD1[2,2:5]) )
lines(data1, col=2, lty=2)
points(data1, pch=22, cex=2, col=2, lty=1)

# set-up labels
labels.text=t(data_cascade)[1,]
labels.col=seq(1,rows_num,1)
# set the axis major tick marks
axisX.at <- 10 ^ c(minX.pow:maxX.pow)
axisY.at <- 10 ^ c(minY.pow:maxY.pow)
# draw the major tick marks and label them using plotmath 
axis(1, at = axisX.at, tcl = -1, labels = parse(text = paste("10^", minX.pow:maxX.pow, sep = "")))
axis(2, at = axisY.at, tcl = -1, labels = parse(text = paste("10^", minY.pow:maxY.pow, sep = "")))
# now do the minor ticks, at 1/10 of each power of 10 interval 
axis(1, at = 1:10 * rep(axisX.at[-1] / 10, each = 10), tcl = -0.5, labels = FALSE)
axis(2, at = 1:10 * rep(axisY.at[-1] / 10, each = 10), tcl = -0.5, labels = FALSE)

lines(t(rbind(c(1.5E5,2.2E5), c(0.5E-3,0.5E-3))), col=1, lty=1)
text(t(c(2.5E5,0.45E-3)), expression(paste(E[a](O[3]) == 1000,K)), adj = c(0,0))
lines(t(rbind(c(1.5E5,2.2E5), c(0.3E-3,0.3E-3))), col=2, lty=2)
text(t(c(2.5E5,0.27E-3)), expression(paste(E[a](O[3]) == 2000,K)), adj = c(0,0), col=2)

points(t(c(6E4,0.17E-3)), pch=8, cex=2, col=1, lty=1)
text(t(c(7.2E4,0.155E-3)), expression(paste(D/H,"=", "0.0001")), adj = c(0,0), col=1, cex=0.9)
points(t(c(6E4,0.1E-3)), pch=5, cex=2, col=1, lty=1)
text(t(c(7.2E4,0.088E-3)), expression(paste(D/H,"=", "0.001")), adj = c(0,0), col=1, cex=0.9)

points(t(c(3E5,0.17E-3)), pch=2, cex=2, col=1, lty=1)
text(t(c(3.6E5,0.155E-3)), expression(paste(D/H,"=", "0.1")), adj = c(0,0), col=1, cex=0.9)
points(t(c(3E5,0.1E-3)), pch=22, cex=2, col=1, lty=1)
text(t(c(3.6E5,0.088E-3)), expression(paste(D/H,"=", "1.0")), adj = c(0,0), col=1, cex=0.9)
